Poverty report
Карина :)
10.12.2022
Предварительный анализ данных
library(ggplot2)
library(dplyr)
library(lubridate)
library(vroom)
library(readxl)
library(emoGG)
library(GGally)
library(psych)
library(psych)
library(nortest)
library(qqplotr)
library(reshape)
library(ppcor)
library(car)
library(lm.beta)
library(ellipse)
library(vegan)
library(olsrr)Считываем данные.
excel_sheets("D:/data/nov/poverty.xls")## [1] "Sheet1"dataf <- read_excel("D:/data/nov/poverty.xls", sheet="Sheet1")Посмотрим на структуру и первые строчки.
head(dataf)## # A tibble: 6 x 9
## ...1 BIRTH DEATH INF_DEAH LIFE_M LIFE_F GNP GROUP COUNTRY
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 ALBANIA 24.7 5.7 30.8 69.6 75.5 600 EAST_E ALBANIA
## 2 BULGARIA 12.5 11.9 14.4 68.3 74.7 2250 EAST_E BULGARIA
## 3 CZECHOSL 13.4 11.7 11.3 71.8 77.7 2980 EAST_E CZECHOSL
## 4 FORMER_E 12 12.4 7.6 69.8 75.9 NA EAST_E FORMER_E
## 5 HUNGARY 11.6 13.4 14.8 65.4 73.8 2780 EAST_E HUNGARY
## 6 POLAND 14.3 10.2 16 67.2 75.7 1690 EAST_E POLANDdim(dataf)## [1] 97 9Это данные по различным показателям уровня жизни для 97 стран. Оставим только один столбец с названием страны, пропуск значения заменен на NA. Группа - это фактор.
datafm <- dataf[ , -9]
names(datafm)[1] <- "COUNTRY"
datafm$GROUP <- as.factor(datafm$GROUP)
head(datafm)## # A tibble: 6 x 8
## COUNTRY BIRTH DEATH INF_DEAH LIFE_M LIFE_F GNP GROUP
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 ALBANIA 24.7 5.7 30.8 69.6 75.5 600 EAST_E
## 2 BULGARIA 12.5 11.9 14.4 68.3 74.7 2250 EAST_E
## 3 CZECHOSL 13.4 11.7 11.3 71.8 77.7 2980 EAST_E
## 4 FORMER_E 12 12.4 7.6 69.8 75.9 NA EAST_E
## 5 HUNGARY 11.6 13.4 14.8 65.4 73.8 2780 EAST_E
## 6 POLAND 14.3 10.2 16 67.2 75.7 1690 EAST_ECOUNTRY - название страны
BIRTH - рождаемость на 1000 человек населения
DEAH - смертность на 1000 населения
INF_DEAH - младенческая смертность в возрасте до 1 года на 1000 населения
LIFE_M - продолжительность жизни мужчин
LIFE_F - продолжительность жизни женщин
GNP - валовой национальный продукт на душу населения в долларах США
GROUP - группа страны (в какой части мира находится)
Построим pairs plot.
ggpairs(datafm, columns = 2:8, diag = list(continuous = "barDiag"))
Есть сильно несимметричные распределения и есть неоднородности (видны несколько облаков точек). Категоризующая переменная, объясняющая эту неоднородность - часть мира, в которой расположена страна. Прологарифмируем GNP и раскрасим по группам. Снова строим pairs plot.
datafm$GNP <- log(datafm$GNP)
ggpairs(datafm, columns = 2:8, aes(color = GROUP, alpha = 0.6), diag = list(continuous = "barDiag"))
Зависимости стали более линейными, а распределения более симметричными.
Разделяем на группы.
data1 <- datafm %>% filter(GROUP == "EAST_E")
data2 <- datafm %>% filter(GROUP == "SOU_A")
data3 <- datafm %>% filter(GROUP == "WEST_E_A")
data4 <- datafm %>% filter(GROUP == "ASIA")
data5 <- datafm %>% filter(GROUP == "MID_EAST")
data6 <- datafm %>% filter(GROUP == "AFRICA")
data7<-dataПосмотрим на descriptive statistics. Она показывает основные характеристики распределений по группам. У этой функции следующие значения:
mean - среднее,
sd - стандартное отклонение,
median - медиана,
trimmed - усеченное среднее,
mad - среднее абсолютное отклонение от медианы,
min - минимум,
max - максимум,
range - размах (разница между максимумом и минимумом),
skew - асиммерия,
kurtosis - эксцесс,
se - стандартная ошибка.
Посмотрим на группу AFRICA, в которой 27 стран.
Коэффициент асимметрии для рождаемости равен -0.94. Значит распределение сильно скошено вправо (также это понятно из того, что среднее меньше медианы). Минимум расположен далеко от среднего, так как это значение маловероятно. Оно находится на расстоянии больше, чем 2 сигма (13.43 > 2*5.69). А вероятность быть на расстоянии от мат.ожидания больше двух сигм равна всего 0.046. Эксцесс равен -0.15, значит пик более пологий, а хвосты тяжелее, чем у нормального распределения.
data_fm <- datafm[ , -1]
data_fm <- data_fm[ , -7]
describeBy(data_fm, datafm$GROUP)##
## Descriptive statistics by group
## group: AFRICA
## vars n mean sd median trimmed mad min max range skew
## BIRTH 1 27 44.53 5.69 46.10 45.03 3.56 31.10 52.20 21.1 -0.94
## DEATH 2 27 14.62 4.80 14.00 14.38 5.49 7.30 25.00 17.7 0.48
## INF_DEAH 3 27 99.79 30.58 103.00 99.83 43.00 49.40 154.00 104.6 0.14
## LIFE_M 4 27 50.64 7.10 50.30 50.57 9.19 38.10 64.90 26.8 0.04
## LIFE_F 5 27 54.14 7.03 53.70 54.26 9.49 41.20 66.40 25.2 -0.10
## GNP 6 27 6.20 1.04 6.04 6.17 0.97 4.38 8.58 4.2 0.28
## kurtosis se
## BIRTH -0.15 1.09
## DEATH -0.85 0.92
## INF_DEAH -1.34 5.89
## LIFE_M -1.00 1.37
## LIFE_F -1.16 1.35
## GNP -0.65 0.20
## ------------------------------------------------------------
## group: ASIA
## vars n mean sd median trimmed mad min max range skew
## BIRTH 1 17 29.62 8.91 30.50 29.97 12.16 11.7 42.20 30.50 -0.25
## DEATH 2 17 10.22 4.66 8.80 10.01 3.85 4.9 18.70 13.80 0.63
## INF_DEAH 3 17 67.72 51.35 64.00 64.24 59.30 6.1 181.60 175.50 0.57
## LIFE_M 4 17 60.49 8.70 62.50 60.87 7.86 41.0 74.30 33.30 -0.60
## LIFE_F 5 17 63.28 10.62 66.10 63.57 9.79 42.0 80.10 38.10 -0.42
## GNP 6 17 6.49 1.35 6.35 6.40 0.72 4.7 9.56 4.86 1.00
## kurtosis se
## BIRTH -1.07 2.16
## DEATH -1.21 1.13
## INF_DEAH -0.88 12.45
## LIFE_M -0.53 2.11
## LIFE_F -1.04 2.58
## GNP 0.15 0.33
## ------------------------------------------------------------
## group: EAST_E
## vars n mean sd median trimmed mad min max range skew kurtosis
## BIRTH 1 11 14.76 3.69 13.60 14.01 1.63 11.6 24.7 13.1 1.69 1.91
## DEATH 2 11 10.55 2.08 10.70 10.78 1.78 5.7 13.4 7.7 -0.86 0.14
## INF_DEAH 3 11 17.37 7.06 14.80 16.97 5.19 7.6 30.8 23.2 0.55 -1.05
## LIFE_M 4 11 67.69 2.14 67.20 67.58 2.08 64.6 71.8 7.2 0.36 -1.08
## LIFE_F 5 11 74.99 1.39 74.80 74.98 1.33 72.4 77.7 5.3 0.06 -0.45
## GNP 6 9 7.48 0.48 7.54 7.48 0.27 6.4 8.0 1.6 -1.03 0.10
## se
## BIRTH 1.11
## DEATH 0.63
## INF_DEAH 2.13
## LIFE_M 0.65
## LIFE_F 0.42
## GNP 0.16
## ------------------------------------------------------------
## group: MID_EAST
## vars n mean sd median trimmed mad min max range skew
## BIRTH 1 11 33.90 8.63 31.70 33.89 13.20 22.30 45.6 23.30 0.00
## DEATH 2 11 6.75 2.65 7.60 6.73 1.78 2.20 11.5 9.30 -0.12
## INF_DEAH 3 11 47.58 30.77 44.00 45.07 40.03 9.70 108.1 98.40 0.43
## LIFE_M 4 11 64.82 5.01 63.10 64.81 2.08 55.80 73.9 18.10 0.19
## LIFE_F 5 11 67.86 6.06 67.00 68.23 3.26 55.00 77.4 22.40 -0.31
## GNP 6 11 8.52 0.90 8.56 8.52 1.09 7.12 9.9 2.77 0.03
## kurtosis se
## BIRTH -1.81 2.60
## DEATH -0.96 0.80
## INF_DEAH -1.08 9.28
## LIFE_M -0.77 1.51
## LIFE_F -0.32 1.83
## GNP -1.39 0.27
## ------------------------------------------------------------
## group: SOU_A
## vars n mean sd median trimmed mad min max range skew
## BIRTH 1 12 29.18 7.39 28.45 28.55 6.60 18.0 46.60 28.60 0.70
## DEATH 2 12 9.42 5.51 7.65 8.54 1.93 4.4 23.20 18.80 1.49
## INF_DEAH 3 12 51.33 31.70 42.50 48.78 29.43 17.1 111.00 93.90 0.81
## LIFE_M 4 12 62.71 4.93 63.40 63.31 3.78 51.0 68.40 17.40 -0.97
## LIFE_F 5 12 68.53 5.31 68.05 69.19 2.97 55.4 75.10 19.70 -0.90
## GNP 6 12 7.26 0.66 7.35 7.34 0.69 5.8 7.89 2.09 -0.78
## kurtosis se
## BIRTH 0.23 2.13
## DEATH 0.85 1.59
## INF_DEAH -0.71 9.15
## LIFE_M 0.15 1.42
## LIFE_F 0.56 1.53
## GNP -0.58 0.19
## ------------------------------------------------------------
## group: WEST_E_A
## vars n mean sd median trimmed mad min max range skew kurtosis
## BIRTH 1 19 12.85 1.94 13.20 12.81 1.93 9.7 16.70 7.00 -0.02 -1.00
## DEATH 2 19 9.43 1.49 9.40 9.45 1.78 6.7 11.90 5.20 -0.08 -1.12
## INF_DEAH 3 19 7.86 1.87 7.50 7.75 0.74 4.5 13.10 8.60 0.95 1.41
## LIFE_M 4 19 71.50 2.66 72.00 71.60 1.93 65.4 75.90 10.50 -0.84 -0.02
## LIFE_F 5 19 78.18 2.30 78.60 78.31 1.93 72.4 81.80 9.40 -0.88 0.24
## GNP 6 18 9.72 0.42 9.82 9.75 0.28 8.7 10.17 1.47 -1.05 -0.09
## se
## BIRTH 0.45
## DEATH 0.34
## INF_DEAH 0.43
## LIFE_M 0.61
## LIFE_F 0.53
## GNP 0.10Минимумы и максимумы, какие конкретно страны.
EAST_E
Рождаемость
Венгрия 11.6 Албания 24.7
Смертность
Албания 5.7 Венгрия 13.4
Смертность младенцев
Неизвестная страна (FORMER_E) 7.6 Албания 30.8
Продолжительность жизни мужчин
СССР 64.6 Чехия 71.8
Продолжительность жизни женщин
Румыния 72.4 Чехия 77.7
ВНП
Албания 1.86 Чехия 2.08
SOU_A
Рождаемость
Уругвай 18 Боливия 46.6
Смертность
Венесуэла 4.4 Мексика 23.2
Смертность младенцев
Чили 17.1 Боливия 111
Продолжительность жизни мужчин
Боливия 51 Уругвай 68.4
Продолжительность жизни женщин
Боливия 55.4 Чили 75.1
ВНП
Албания 1.76 Чехия 2.07
WEST_E_A
Рождаемость
Италия 9.7 США 16.7
Смертность
Япония 6.7 Дания 11.9
Смертность младенцев
Япония 4.5 Португалия 13.1
Продолжительность жизни мужчин
Греция 64.4 Япония 75.9
Продолжительность жизни женщин
Португалия 72.4 Япония 81.8
ВНП
Греция 2.16 Финляндия 2.32
ASIA
Рождаемость
Гонконг 11.7 Бангладеш 42.2
Смертность
Гонконг 4.9 Афганистан 18.7
Смертность младенцев
Гонконг 6.1 Афганистан 182
Продолжительность жизни мужчин
Афганистан 41 Гонконг 74.3
Продолжительность жизни женщин
Афганистан 42 Гонконг 80.1
ВНП
Монголия 1.55 Гонконг 2.26
MID_EAST
Рождаемость
Израиль 22.3 Оман 45.6
Смертность
Кувейт 2.2 Иран 11.5
Смертность младенцев
Израиль 9.7 Иран 108
Продолжительность жизни мужчин
Иран 55.8 Израиль 73.9
Продолжительность жизни женщин
Афганистан 55 Гонконг 77.4
ВНП
Иордания 1.96 Неизвестная страна (UNITED_A) 2.29
AFRICA
Рождаемость
Тунис 31.1 Гана 44.4
Смертность
Тунис 7.3 Намибия 12.1
Смертность младенцев
Египет 49.4 Ливия 82
Продолжительность жизни мужчин
Малави 38.1 Нигерия 48.8
Продолжительность жизни женщин
Малави 41.2Нигерия 52.2
ВНП
Мозамбик 1.48 Нигерия 1.77
О виде распределений и о сравнении распределений
Выясним, близки ли распределения к нормальным. Будем анализировать вид распределений признаков по группам, разделяем по частям мира.
Normal probability plot.
df_new <- datafm %>%
as.data.frame(datafm) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df1_new <- data1 %>%
as.data.frame(data1) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df2_new <- data2 %>%
as.data.frame(data2) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df3_new <- data3 %>%
as.data.frame(data3) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df4_new <- data4 %>%
as.data.frame(data4) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df5_new <- data5 %>%
as.data.frame(data5) %>%
melt(id.vars = c("COUNTRY", "GROUP"))
df6_new <- data6 %>%
as.data.frame(data6) %>%
melt(id.vars = c("COUNTRY", "GROUP"))Отклонения от прямой линии указывают на отклонения от нормальности.
Далее будем сравнивать распределения признаков у SOU_A и MID_EAST, поэтому построим для них normal probability plot.
Видим, что у распределений Рождаемости, Смертности младенцев для группы MID_EAST хвосты более легкие, чем у нормального распределения. (Отличием от тяжелых является направление отклонения от прямой линии для нескольких первых и нескольких последних точек. Для легких хвостов первые несколько точек показывают отклонение от линии над линией, а последние несколько точек показывают отклонение от прямой линии под линией. Для тяжелых хвостов эта схема обратная.)
Normal probability plot для MID_EAST.
ggplot(df5_new, aes(sample = value)) +
stat_qq_point(size = 2) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_qq_line(color = "darkorchid4")+
ggtitle("Normal probability plot MID_EAST")
Normal probability plot для SOU_A.
ggplot(df2_new, aes(sample = value)) +
stat_qq_point(size = 2) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_qq_line(color = "darkorchid4")+
ggtitle("Normal probability plot SOU_A")
Распределения похожи на нормальные.
PP-plot.
То есть опять используем графический метод для оценки того, следует ли набор данных заданному распределению (нормальному).
Отклонения от прямой линии указывают на отклонения от нормального распределения.
PP-plot для SOU_A.
ggplot(df2_new, aes(sample = value)) +
stat_pp_point(size = 2) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_pp_line(color = "blue")+
ggtitle("PP-plot MID_EAST")
PP-plot для MID_EAST.
ggplot(df5_new, aes(sample = value)) +
stat_pp_point(size = 2) +
facet_wrap(~variable, scales = "free") +
labs(x = "", y = "") +
stat_pp_line(color = "blue")+
ggtitle("PP-plot MID_EAST")
Распределения похожи на нормальные.
Проверка на нормальность по критериям.
Будем рассматривать отдельно по группам, в какой части мира расположена страна. Проверяем по критериям Лиллиефорса (модификация критерия Колмогорова-Смирнова), Андерсона-Даринга (критерий типа омега-квадрат, придает хвостам больший вес, чем тест KS), критерий хи-квадрат Пирсона для сложной гипотезы нормальности, Шапиро-Уилка (примерно квадрат корреляции между x и y в normal probability ploy).
Не будем проводить тесты для всех признаков всех групп. Рассматривается по 2 группы для признаков Рождаемость, Смертность младенцев, Продолжительность жизни мужчин.
Рождаемость
WEST_E_A
lillie.test(data3$BIRTH)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data3$BIRTH
## D = 0.097305, p-value = 0.9069ad.test(data3$BIRTH)##
## Anderson-Darling normality test
##
## data: data3$BIRTH
## A = 0.20873, p-value = 0.84pearson.test(data3$BIRTH)##
## Pearson chi-square normality test
##
## data: data3$BIRTH
## P = 3.4737, p-value = 0.4819shapiro.test(data3$BIRTH)##
## Shapiro-Wilk normality test
##
## data: data3$BIRTH
## W = 0.97024, p-value = 0.7813Гипотеза не отвергается при стандартных уровнях значимости.
SOU_A
lillie.test(data2$BIRTH)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data2$BIRTH
## D = 0.17612, p-value = 0.3835ad.test(data2$BIRTH)##
## Anderson-Darling normality test
##
## data: data2$BIRTH
## A = 0.40813, p-value = 0.2916pearson.test(data2$BIRTH)##
## Pearson chi-square normality test
##
## data: data2$BIRTH
## P = 8, p-value = 0.04601shapiro.test(data2$BIRTH)##
## Shapiro-Wilk normality test
##
## data: data2$BIRTH
## W = 0.92767, p-value = 0.356При уровне значимостии 0.05 гипотеза отвергается по критерию хи-квадрат и не отвергается по остальным критериям.
Смертность младенцев
SOU_A
lillie.test(data2$INF_DEAH)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data2$INF_DEAH
## D = 0.18964, p-value = 0.273ad.test(data2$INF_DEAH)##
## Anderson-Darling normality test
##
## data: data2$INF_DEAH
## A = 0.65002, p-value = 0.06734pearson.test(data2$INF_DEAH)##
## Pearson chi-square normality test
##
## data: data2$INF_DEAH
## P = 4, p-value = 0.2615shapiro.test(data2$INF_DEAH)##
## Shapiro-Wilk normality test
##
## data: data2$INF_DEAH
## W = 0.8588, p-value = 0.04722При уровне значимостии 0.05 гипотеза отвергается по критерию Шапиро-Уилка и не отвергается по остальным критериям.
MID_EAST
lillie.test(data5$INF_DEAH)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data5$INF_DEAH
## D = 0.13094, p-value = 0.8625ad.test(data5$INF_DEAH)##
## Anderson-Darling normality test
##
## data: data5$INF_DEAH
## A = 0.28767, p-value = 0.5492pearson.test(data5$INF_DEAH)##
## Pearson chi-square normality test
##
## data: data5$INF_DEAH
## P = 2.6364, p-value = 0.4512shapiro.test(data5$INF_DEAH)##
## Shapiro-Wilk normality test
##
## data: data5$INF_DEAH
## W = 0.93873, p-value = 0.5057Гипотеза не отвергается при стандартных уровнях значимости.
Продолжительность жизни мужчин
SOU_A
lillie.test(data2$LIFE_M)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data2$LIFE_M
## D = 0.20085, p-value = 0.1995ad.test(data2$LIFE_M)##
## Anderson-Darling normality test
##
## data: data2$LIFE_M
## A = 0.45735, p-value = 0.217pearson.test(data2$LIFE_M)##
## Pearson chi-square normality test
##
## data: data2$LIFE_M
## P = 1, p-value = 0.8013shapiro.test(data2$LIFE_M)##
## Shapiro-Wilk normality test
##
## data: data2$LIFE_M
## W = 0.89971, p-value = 0.1572Гипотеза не отвергается при стандартных уровнях значимости.
MID_EAST
lillie.test(data5$LIFE_M)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: data5$LIFE_M
## D = 0.18543, p-value = 0.3617ad.test(data5$LIFE_M)##
## Anderson-Darling normality test
##
## data: data5$LIFE_M
## A = 0.38357, p-value = 0.3309pearson.test(data5$LIFE_M)##
## Pearson chi-square normality test
##
## data: data5$LIFE_M
## P = 2.6364, p-value = 0.4512shapiro.test(data5$LIFE_M)##
## Shapiro-Wilk normality test
##
## data: data5$LIFE_M
## W = 0.94745, p-value = 0.6116Гипотеза не отвергается при стандартных уровнях значимости.
Как по результатам проверки на нормальность прикинуть, значения асимметрии и эксцесса были около нуля или существенными? Если гипотеза о нормальности отвергается, то значения асимметрии и эксцесса можно считать существенными.
Ящики с усами.
Сравненим распределения в группах с помощью ящиков с усами.
cols <- c("#7edecc", "#79c4db", "#aede81", "#ded17e", "#8d76e8", "#de7eb6")
ggplot(datafm, aes(x = GROUP, y = BIRTH, fill = GROUP)) +
stat_boxplot(geom = "errorbar",
width = 0.25) +
geom_boxplot(alpha = 0.8,
colour = "#474747",
outlier.colour = 1) +
scale_fill_manual(values = cols) +
ggtitle("Рождаемость")
ggplot(datafm, aes(x = GROUP, y = DEATH, fill = GROUP)) +
stat_boxplot(geom = "errorbar",
width = 0.25) +
geom_boxplot(alpha = 0.8,
colour = "#474747",
outlier.colour = 1) +
scale_fill_manual(values = cols) +
ggtitle("Смертность")
ggplot(datafm, aes(x = GROUP, y = INF_DEAH, fill = GROUP)) +
stat_boxplot(geom = "errorbar",
width = 0.25) +
geom_boxplot(alpha = 0.8,
colour = "#474747",
outlier.colour = 1) +
scale_fill_manual(values = cols) +
ggtitle("Смертность младенцев")
ggplot(datafm, aes(x = GROUP, y = LIFE_M, fill = GROUP)) +
stat_boxplot(geom = "errorbar",
width = 0.25) +
geom_boxplot(alpha = 0.8,
colour = "#474747",
outlier.colour = 1) +
scale_fill_manual(values = cols) +
ggtitle("Продолжительность жизни мужчин")
ggplot(datafm, aes(x = GROUP, y = LIFE_F, fill = GROUP)) +
stat_boxplot(geom = "errorbar",
width = 0.25) +
geom_boxplot(alpha = 0.8,
colour = "#474747",
outlier.colour = 1) +
scale_fill_manual(values = cols) +
ggtitle("Продолжительность жизни женщин")
t-критерий.
Смертность младенцев.
Сначала проверим гипотезу о равенстве дисперсий для двух распределений. Рассматриваем Смертность младенцев для SOU_A и MID_EAST.
Критерий Фишера можно использовать только для нормальных распределений, а один из тестов отверг гипотезу о нормальности Смертности младенцев в SOU_A, поэтому будем использовать критерий Левена.
data25 <- datafm %>% filter(GROUP == "SOU_A" | GROUP == "MID_EAST" )
leveneTest(INF_DEAH ~ GROUP, data25, center = "mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 1 0 0.998
## 21Гипотеза не отвергается при стандартных уровнях значимости. Но нельзя считать, что дисперсии равны.
Будем использовать двухвыборочный t-критерий для независимых выборок с равными дисперсиями и с различными дисперсиями, критерий асимптотический. t-критерий точный, когда данные нормальные.
Если бы был сбалансированный дизайн, то не важно равны дисперсии или нет. (Но он не сбалансирован, разные объемы выборки).
t.test(data2$INF_DEAH, data5$INF_DEAH, var.equal=TRUE)##
## Two Sample t-test
##
## data: data2$INF_DEAH and data5$INF_DEAH
## t = 0.28688, df = 21, p-value = 0.777
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -23.39176 30.87813
## sample estimates:
## mean of x mean of y
## 51.32500 47.58182t.test(data2$INF_DEAH, data5$INF_DEAH, var.equal=FALSE)##
## Welch Two Sample t-test
##
## data: data2$INF_DEAH and data5$INF_DEAH
## t = 0.28726, df = 20.921, p-value = 0.7767
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -23.36151 30.84787
## sample estimates:
## mean of x mean of y
## 51.32500 47.58182Гипотезы о равенстве математических ожиданий не отвергается при стандартных уровнях значимости.
Продолжительность жизни мужчин.
Аналогично посмотрим на другой признак - Продолжительность жизни мужчин.
Сначала проверяем гипотезу о равенсте дисперсий. Здесь можно использовать критерий Фишера.
var.test(data2$LIFE_M, data5$LIFE_M)##
## F test to compare two variances
##
## data: data2$LIFE_M and data5$LIFE_M
## F = 0.9652, num df = 11, denom df = 10, p-value = 0.9474
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.2633614 3.4029667
## sample estimates:
## ratio of variances
## 0.9651967Гипотеза не отвергается при стандартных уровнях значимости.
Будем использовать двухвыборочный t-критерий для независимых выборок с равными и различными дисперсиями.
t.test(data2$LIFE_M, data5$LIFE_M, var.equal=TRUE)##
## Two Sample t-test
##
## data: data2$LIFE_M and data5$LIFE_M
## t = -1.0175, df = 21, p-value = 0.3205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.422153 2.202456
## sample estimates:
## mean of x mean of y
## 62.70833 64.81818t.test(data2$LIFE_M, data5$LIFE_M, var.equal=FALSE)##
## Welch Two Sample t-test
##
## data: data2$LIFE_M and data5$LIFE_M
## t = -1.0167, df = 20.754, p-value = 0.321
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.428751 2.209054
## sample estimates:
## mean of x mean of y
## 62.70833 64.81818Гипотезы о равенстве математических ожиданий не отвергается при стандартных уровнях значимости.
Критерий Вилкоксона.
Критерий Вилкоксона непараметрический. Хорош тем, что можно использовать при малых объемах выборки и робастностью (свойство статистического метода, характеризующее независимость влияния выбросов на результат исследования).
Плох тем, что у него мощность меньше, чем у t-критерия.
У нас очень маленькие объемы выборки. Для SOU_A - это 12 индивидов, для MID_EAST - 11 индивидов. Используем критерий Вилкоксона.
wilcox.test(data2$INF_DEAH, data5$INF_DEAH)##
## Wilcoxon rank sum test with continuity correction
##
## data: data2$INF_DEAH and data5$INF_DEAH
## W = 69.5, p-value = 0.8534
## alternative hypothesis: true location shift is not equal to 0wilcox.test(data2$LIFE_M, data5$LIFE_M)##
## Wilcoxon rank sum test with continuity correction
##
## data: data2$LIFE_M and data5$LIFE_M
## W = 56, p-value = 0.5587
## alternative hypothesis: true location shift is not equal to 0Гипотезы не отвергаются при стандартных уровнях значимости.
Гипотеза о том, что P(кси1 > кси2) = P(кси1 < кси2). t-критерий, но примененный к рангам.
Критерий Колмогорова Смирнова.
С помощью критерия Колмогорова-Смирнова можно сравнивать распределения в целом (умеет сравнивать формы распределений).
ks.test(data2$INF_DEAH, data5$INF_DEAH)##
## Two-sample Kolmogorov-Smirnov test
##
## data: data2$INF_DEAH and data5$INF_DEAH
## D = 0.27273, p-value = 0.7868
## alternative hypothesis: two-sidedks.test(data2$LIFE_M, data5$LIFE_M)##
## Two-sample Kolmogorov-Smirnov test
##
## data: data2$LIFE_M and data5$LIFE_M
## D = 0.27273, p-value = 0.7868
## alternative hypothesis: two-sidedГипотезы о равенстве распределений не отвергаются при стандартных уровнях значимости.
Об анализе зависимостей
У нас неоднородные данные, будем изучать зависимости только внутри групп по-отдельности. Построим сначала pairs plot для группы AFRICA.
Коэффициент корреляции Спирмана не реагирует на монотонные преобразования и почти не реагирует на выбросы. Посмотрим на смертность и продолжительность жизни мужчин, видим, что есть выброс. Значит, коэффициент корреляции Спирмана должен быть больше, чем коэффициент корреляции Пирсона. Проверим это, построив соответствующие матрицы.
ggpairs(data6, columns = 2:8, diag = list(continuous = "barDiag"))
Матрица корреляций Пирсона
dplyr::select(data6, -COUNTRY, -GROUP) %>%
cor(method = "pearson", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 3))) +
scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) +
ggtitle("Pearson AFRICA") +
theme(axis.text.x = element_text(angle = 50, hjust = 1))
Коэффициент корреляции Пирсона для смертности и продолжительности жизни мужчин равен -0.935.
Матрица корреляций Спирмана
dplyr::select(data6, -COUNTRY, -GROUP) %>%
cor(method = "spearman", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
geom_raster(aes(fill = value)) +
geom_text(aes(label = round(value, 3))) +
scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) +
ggtitle("Spearman AFRICA") +
theme(axis.text.x = element_text(angle = 50, hjust = 1))
Коэффициент корреляции Спирмана для смертности и продолжительности жизни мужчин равен -0.951.
Для нормальных распределений коэффициенты Пирсона и Спирмана примерно равны.
Видим достаточно большую корреляцию между продолжительностью жизни (и мужчин, и женщин) и ВНП. Что является причиной, а что следствием? Наверно, чем больше ВНП страны на душу населения, тем больше будет продолжительности жизни людей в этой стране.
Посмотрим на частную корреляцию между Рождаемостью и Смертностью за вычетом ВНП.
(data6 %>%
dplyr::select(BIRTH, DEATH, GNP) %>%
pcor(method = "pearson"))$estimate["BIRTH", "DEATH"]## [1] 0.3536362Частная корреляция 0.3536 между ними меньше, чем обычная корреляция 0.609.
Посмотрим на частную корреляцию между Продолжительностью жизни мужчин и ВНП за вычетом Смертности.
(data6 %>%
dplyr::select(LIFE_M, GNP, DEATH) %>%
pcor(method = "pearson"))$estimate["LIFE_M", "GNP"]## [1] 0.02495314Частная корреляция 0.02495 между ними сильно меньше, чем обычная корреляция 0.658.
Линейная регрессия
#data1 <- datafm %>% filter(GROUP == "EAST_E")
data_asia_mid_east <- datafm %>% filter((GROUP == "ASIA")|(GROUP == "MID_EAST"))
data_africa <- datafm %>% filter(GROUP == "AFRICA")
#head(data_asia_mid_east)ggpairs(data_asia_mid_east, columns = 2:8, diag = list(continuous = "barDiag"))
ggpairs(data_africa, columns = 2:8, diag = list(continuous = "barDiag"))
model_1 <- lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa)
summary(lm.beta(model_1))##
## Call:
## lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1095 -1.9535 -0.2117 1.8078 8.2824
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 76.04748 NA 10.45968 7.271 2.12e-07 ***
## BIRTH -0.39253 -0.31450 0.14936 -2.628 0.015 *
## INF_DEAH -0.13642 -0.58793 0.02794 -4.883 6.24e-05 ***
## GNP 0.91661 0.13486 0.85326 1.074 0.294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.398 on 23 degrees of freedom
## Multiple R-squared: 0.7972, Adjusted R-squared: 0.7708
## F-statistic: 30.14 on 3 and 23 DF, p-value: 3.821e-08ellipse_1 <- ellipse(model_1, which=c(2,3), npoints=1000)
#ellipse_1
ggplot() + geom_point(aes(x = ellipse_1[,1], y = ellipse_1[,2]))+xlab("BIRTH")+ylab("INF_DEAH")
ggcorr(data_africa,
label = T,
label_size = 4,
label_round = 2,
hjust = 1,
size = 5,
color = "royalblue",
layout.exp = 5,
low = "green3",
mid = "gray95",
high = "darkorange",
name = "Correlation")
ols_step_backward_p(model_1)## [1] "No variables have been removed from the model."#step(model_1, direction = "backward")ols_step_forward_p(model_1)##
## Selection Summary
## -------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------
## 1 INF_DEAH 0.6903 0.6779 12.1276 155.7736 4.0274
## 2 BIRTH 0.7871 0.7693 3.1540 147.6611 3.4085
## 3 GNP 0.7972 0.7708 4.0000 148.3393 3.3976
## -------------------------------------------------------------------------#step(model_1, direction = "forward")shapiro.test(model_1$residuals)##
## Shapiro-Wilk normality test
##
## data: model_1$residuals
## W = 0.98271, p-value = 0.9178ggplot(data_africa ,aes(x = model_1$fitted.values, y = model_1$residuals)) +
geom_point()+
xlab("Predicted") +
ylab("Residuals")
data_africaa <- data_africa[ , -1]
data_africaa <- data_africaa[ , -7]
datt <- data.frame('n' = 1:27, 'country' = data_africa$COUNTRY,
'mahalanobis' = mahalanobis(data_africaa,
colMeans(data_africaa),
cov(data_africaa)))
ggplot(datt, aes(x = n, y = mahalanobis)) +
geom_point(size= 2) +
xlab("n") +
ylab("Mahalanobis' distance")
dat <- data.frame('n' = 1:27, 'country' = data_africa$COUNTRY,
'cook' = cooks.distance(model_1))
ggplot(dat, aes(x = n, y = cook)) +
geom_point(size= 2) +
xlab("n") +
ylab("Cook's distance")
ols_plot_resid_stud(model_1)
ols_plot_resid_stud_fit(model_1)
ggplot(data_frame(residuals=rstandard(model_1), studres=studres(model_1)), aes(x=residuals, y=studres))+geom_point()+ geom_abline(scope=1, intercept=0, color = "red")
data_africa_new <- data_africa[ -c(11,15,19),]
model_2 <- lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa_new)
summary(lm.beta(model_2))##
## Call:
## lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa_new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5260 -1.0234 0.1219 1.1299 3.9726
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 95.50820 NA 9.48892 10.065 2.84e-09 ***
## BIRTH -0.49199 -0.37688 0.12434 -3.957 0.000778 ***
## INF_DEAH -0.17736 -0.77363 0.02377 -7.462 3.36e-07 ***
## GNP -0.95544 -0.12095 0.85150 -1.122 0.275125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.578 on 20 degrees of freedom
## Multiple R-squared: 0.885, Adjusted R-squared: 0.8677
## F-statistic: 51.29 on 3 and 20 DF, p-value: 1.421e-09#new_africa <- data.frame("COUNTRY" = "Land", BIRTH = 38, DEATH = 20, #INF_DEAH = 137, LIFE_M = 60, LIFE_F=70, GNP = 700, GROUP="AFRICA")
new_africa2 <- data.frame(COUNTRY = "Land", BIRTH = 38, INF_DEAH = 137, GNP = 700)
new_africa2## COUNTRY BIRTH INF_DEAH GNP
## 1 Land 38 137 700predict(model_1, newdata = new_africa2, interval="confidence")## fit lwr upr
## 1 684.0709 -540.6429 1908.785predict(model_1, newdata = new_africa2, interval="prediction")## fit lwr upr
## 1 684.0709 -540.663 1908.805